
clear
clear all
clear mata
clear matrix
set more off
set matsize 11000
set maxvar 30000
cap log off
capture log close
set emptycells drop
pause on


***********************************************
* USER DEFINE FILEPATH FOR REPLICATION FOLDER *
***********************************************

global path "C:\Users\wb520443\Dropbox\Research Projects\Global pollution\2_analysis\Replication"



***********************************************


global opts		a f plain coll(none) nodep nomti c(b(star fmt(%9.3f)) se(abs par fmt(%9.3f))  ) star(* .10 ** .05 *** .01) noobs nocons

global opts_fs		a f plain coll(none) nodep nomti c(b(star fmt(%9.3f)) se(abs par fmt(%9.3f))) star(* .10 ** .05 *** .01) noobs nocons


global input "$path\data"
global figdat "$path\figdat"
global figures "$path\figures"
global tables "$path\tables"



cap mat define countryfigure=J(105,10,.)
		local row=1
		
	use "$input/global_pm_income.dta"
		replace ISO="IND/PAK" if ISO=="IND"
		rename ISO country
			
		save "$input/temp_income.dta", replace
	
		clear	
	
	use "$input/rwi_pollution_analysis_collapsed.dta", replace
	
	replace error=1/error
	g both=population*error
	g both2=country_population*error
	
	
	cap encode country, g(country_code)
		label save using "$input/labels.do", replace

		
	//Mean, mean_ss, and max refer to how PM25 across years and months is collapse (e.g. max is the max pollution level for a RWI point across all years and months in the sample).
	forvalues i=1(1)103{
		
		cap reghdfe mean_pm25 rwi [pweight=both] if country_code==`i', noabsorb
			cap mat def countryfigure[`row',1]=_b[rwi]
			cap mat def countryfigure[`row',2]=_se[rwi]
		cap reghdfe mean_pm25_ss rwi  [pweight=both]if country_code==`i', noabsorb
			cap mat def countryfigure[`row',3]=_b[rwi]
			cap mat def countryfigure[`row',4]=_se[rwi]
		
		local row=`row'+1
	}	

clear

	
	
clear
			
		svmat2 countryfigure

		drop if countryfigure1==.
		
		g country_code=_n
			run "$input/labels.do"
		
		g HI95_mean=(1.96*countryfigure2)+countryfigure1
		g LI95_mean=countryfigure1-(1.96*countryfigure2)

		g HI99_mean=(2.06*countryfigure2)+countryfigure1
		g LI99_mean=countryfigure1-(2.06*countryfigure2)
		
		
		g HI95_mean_ss=(1.96*countryfigure4)+countryfigure3
		g LI95_mean_ss=countryfigure3-(1.96*countryfigure4)

		g HI99_mean_ss=(2.06*countryfigure4)+countryfigure3
		g LI99_mean_ss=countryfigure3-(2.06*countryfigure4)
		
				
		set scheme plotplain
		
		
		rename (countryfigure1 countryfigure3 countryfigure5 countryfigure7 countryfigure9 ) (mean_beta mean_ss_beta max_beta urban_beta rural_beta)
		

		
		

		label val country_code country_code
			decode country_code, g(country)
			replace country="IND/PAK" if country=="IND_PAK"
			
			
		merge 1:1 country using "$input/temp_income.dta", nogen keep(1 3)
		
		g bin=1 if mean_gdp<2000 
		replace bin=2 if mean_gdp<4000 & mean_gdp>=2000
		replace bin=3 if mean_gdp<6000 & mean_gdp>=4000
		replace bin=4 if mean_gdp<8000 & mean_gdp>=6000
		replace bin=5 if mean_gdp>=8000
		
		gsort - bin - mean_beta
		
		g n=_n
		g n2=n+0.1
		
		

		drop if mean_gdp==.
			labmask n, values(country)
		
		preserve
			drop if mean_beta<-5

			g missing=.
			twoway 	(scatter  n missing , yaxis(1)) ///
						(rspike HI99_mean LI99_mean n , lc(gs10) horizontal yaxis(2)) ///
						(scatter n mean_beta  if mean_gdp<2000, mc("178 34 34")  ms(D)  msize(large) yaxis(2)) ///
						(scatter n mean_beta  if mean_gdp<4000 & mean_gdp>=2000, mc("197 80 38") ms(D) msize(large) yaxis(2)) ///
						(scatter n mean_beta  if mean_gdp<6000 & mean_gdp>=4000, mc("215 118 39") ms(D) msize(large) yaxis(2)) ///
						(scatter n mean_beta  if mean_gdp<8000 & mean_gdp>=6000, mc("232 155 36") ms(D) msize(large) yaxis(2)) ///
						(scatter n mean_beta  if mean_gdp>=8000, mc("255 215 0") ms(D) msize(large) yaxis(2) ///
						ytitle("",axis(1)) yla(none,nogrid axis(1)) ysc(noline axis(1)) ysc(noline axis(2))  ylabel("",alternate angle(0) labsize(small) nogrid axis(2))  ytitle("") xtitle("Correlation between wealth" "and average pollution" " ",size(vlarge)) ysize(9) xsize(2) legend(col(1)  ring(0) position(2) order(3 4 5 6 7) label(3 "<$2000") label(4 "$2000-3999") label(5 "$4000-5999") label(6 "$6000-7999") label(7 ">$8000") size(*2.25)) xlabel(,labsize(large)) xline(0))
						
						
			
						
						graph export "$figures/country_figure_mean_collapsed_income.png", as(png) replace
		restore
		
		drop n
		gsort - bin - mean_ss_beta
		
		g n=_n			
		labmask n, values(country)
		drop if mean_ss_beta>10
		twoway 	(rspike HI99_mean_ss LI99_mean_ss n , lc(gs10) horizontal) ///
					(scatter n mean_ss_beta  if mean_gdp<2000, mc("178 34 34")  ms(D ) msize(large)) ///
					(scatter n mean_ss_beta  if mean_gdp<4000 & mean_gdp>=2000, mc("197 80 38") ms(D) msize(large)) ///
					(scatter n mean_ss_beta  if mean_gdp<6000 & mean_gdp>=4000, mc("215 118 39") ms(D) msize(large)) ///
					(scatter n mean_ss_beta  if mean_gdp<8000 & mean_gdp>=6000, mc("232 155 36") ms(D) msize(large)) ///
					(scatter n mean_ss_beta  if mean_gdp>=8000, mc("255 215 0") ms(D) msize(large) ///
					ytitle("",axis(1)) yla(none,nogrid axis(1)) ysc(noline axis(1))   ytitle("") xtitle("Correlation between wealth" "and average pollution" "with dust and sea salt removed",size(vlarge)) ysize(9) xsize(2) legend(off) xlabel(,labsize(large)) xline(0))	
					
					
					
					
					graph export "$figures/country_figure_mean_collapsed_income_ss.png", as(png) replace
clear
	
	

cap mat define countryfigure=J(105,10,.)
		local row=1
		
	import delimited "$input/country_ID_crosswalk_for_EDGAR.csv"
		rename x ID_0
		g id=_n
	save "$input/edgar_country_crosswalk.dta", replace
	clear

foreach e in AGS  AWB  CHE ENE FFF  IND IRO MNM NFE NMM PRO  RCO REF_TRF SWD_LDF SWD_INC   TNR_Aviation_CDS TNR_Aviation_LTO TNR_Aviation_CRS TNR_Other TNR_Ship TRO_noRES TRO_RES{
			import delimited "$input/all_countries_edgar_`e'.csv", varnames(1)
			rename * `e'
			g id=_n
			
			cap destring `e', replace force
			
			save "$input/all_countries_edgar_`e'.dta", replace
			clear
		}
		
		use "$input/all_countries_edgar_AGS.dta"
		foreach e in  AWB  CHE ENE FFF  IND IRO MNM NFE NMM PRO  RCO REF_TRF SWD_LDF SWD_INC    TNR_Aviation_CDS TNR_Aviation_LTO TNR_Aviation_CRS TNR_Other TNR_Ship TRO_noRES TRO_RES{
			merge 1:1 id using "$input/all_countries_edgar_`e'.dta", nogen keep(1 3)
			
			erase "$input/all_countries_edgar_`e'.dta"
		}
		
		merge 1:1 id using "$input/edgar_country_crosswalk.dta", nogen keep(1 3)
		
		g power=ENE+PRO+RCO
		g industry=REF_TRF+IND+CHE+IRO+NFE
		g transport=TNR_Aviation_CDS +TNR_Aviation_LTO+ TNR_Aviation_CRS +TNR_Other+ TNR_Ship +TRO_noRES+TRO_RES
		g other=SWD_INC+SWD_LDF+FFF
		g agri=MNM+AWB+AGS
		
		keep ID_0 power industry transport other agri
		
		replace ID_0=105 if ID_0==171
		
		g country_code=ID_0
		
		foreach var in power industry transport other agri{
			bysort country_code: egen total_`var'=total(`var')
		}
	
		egen max_source=rowmax(total_power total_industry total_transport total_other total_agri)
		g main_source="power" if max_source==total_power	
			replace main_source="industry" if max_source==total_industry
			replace main_source="transport" if max_source==total_transport
			replace main_source="other" if max_source==total_other
			replace main_source="agri" if max_source==total_agri
		
		duplicates drop ID_0, force
		
		merge 1:1 ID_0 using "$input/iso_id_0_crosswalk.dta"	, nogen keep(1 3)
		
		replace ISO="IND/PAK" if ISO=="IND"
		rename ISO country
			
		save "$input/all_countries_edgar_total.dta", replace
		
		clear	
	
	use "$input/rwi_pollution_analysis_collapsed.dta", replace
	
	replace error=1/error
	g both=population*error
	g both2=country_population*error
	
	
	cap encode country, g(country_code)
		label save using "$input/labels.do", replace
	
		
	//Mean, mean_ss, and max refer to how PM25 across years and months is collapse (e.g. max is the max pollution level for a RWI point across all years and months in the sample).
	forvalues i=1(1)103{
		
		cap reghdfe mean_pm25 rwi [pweight=both] if country_code==`i', noabsorb
			cap mat def countryfigure[`row',1]=_b[rwi]
			cap mat def countryfigure[`row',2]=_se[rwi]
		cap reghdfe mean_pm25_ss rwi  [pweight=both]if country_code==`i', noabsorb
			cap mat def countryfigure[`row',3]=_b[rwi]
			cap mat def countryfigure[`row',4]=_se[rwi]
		
		local row=`row'+1
	}	
clear
	
	
	
clear
			
		svmat2 countryfigure

		drop if countryfigure1==.
		
		g country_code=_n
			run "$input/labels.do"
		
		g HI95_mean=(1.96*countryfigure2)+countryfigure1
		g LI95_mean=countryfigure1-(1.96*countryfigure2)

		g HI99_mean=(2.06*countryfigure2)+countryfigure1
		g LI99_mean=countryfigure1-(2.06*countryfigure2)
		
		
		g HI95_mean_ss=(1.96*countryfigure4)+countryfigure3
		g LI95_mean_ss=countryfigure3-(1.96*countryfigure4)

		g HI99_mean_ss=(2.06*countryfigure4)+countryfigure3
		g LI99_mean_ss=countryfigure3-(2.06*countryfigure4)
		
				
		set scheme plotplain
		
		
		rename (countryfigure1 countryfigure3 countryfigure5 countryfigure7 countryfigure9 ) (mean_beta mean_ss_beta max_beta urban_beta rural_beta)
		
		label val country_code country_code
			decode country_code, g(country)
			replace country="IND/PAK" if country=="IND_PAK"
		
		merge 1:1 country using "$input/all_countries_edgar_total.dta", nogen keep(1 3)
		
		g bin=1 if main_source=="power"
		replace bin=2 if main_source=="industry"
		replace bin=3 if main_source=="transport"
		replace bin=4 if main_source=="agri"
		replace bin=5 if main_source=="other"
		
		gsort - bin - mean_beta
		g n=_n
		g n2=n+0.1
		
		

		
			labmask n, values(country)
			
		
		preserve
		drop if mean_beta<-5
			g missing=.
			twoway 	(scatter  n missing , yaxis(1)) /// 
						(rspike HI99_mean LI99_mean n , lc(gs10) horizontal yaxis(2)) ///
						(scatter n mean_beta  if main_source=="power", mc("0 0 128")  ms(D) msize(large) yaxis(2)) ///
						(scatter n mean_beta  if main_source=="industry", mc("0 54 89") ms(D) msize(large) yaxis(2)) ///
						(scatter n mean_beta  if main_source=="transport", mc("26 77 82") ms(D) msize(large) yaxis(2)) ///
						(scatter n mean_beta  if main_source=="agri", mc("34 101 76") ms(D) msize(large) yaxis(2)) ///
						(scatter n mean_beta  if main_source=="other", mc("33 126 55") ms(D) msize(large) yaxis(2) ///
						ytitle("",axis(1)) yla(none,nogrid axis(1)) ysc(noline axis(1)) ysc(noline axis(2))  ylabel("",alternate angle(0) labsize(small) nogrid axis(2))  ytitle("") xtitle("Correlation between wealth" "and average pollution" " ", size(vlarge)) ysize(9) xsize(2) legend(col(1)  ring(0) position(2) order(3 4 5 6 ) label(3 "Power sector") label(4 "Industry") label(5 "Transportation") label(6 "Agriculture") size(*2.25)) xlabel(,labsize(large)) xline(0))
						
						
						graph export "$figures/country_figure_mean_collapsed_edgar.png", as(png) replace
		
		restore			
		
		drop n			
		gsort - bin - mean_ss_beta
		g n=_n	
		labmask n, values(country)
		drop if mean_ss_beta>10
		twoway 	(rspike HI99_mean_ss LI99_mean_ss n , lc(gs10) horizontal) ///
					(scatter n mean_ss_beta  if main_source=="power", mc("0 0 128")  ms(D) msize(large)) ///
					(scatter n mean_ss_beta  if main_source=="industry", mc("0 54 89") ms(D) msize(large)) ///
					(scatter n mean_ss_beta  if main_source=="transport", mc("26 77 82") ms(D) msize(large)) ///
					(scatter n mean_ss_beta  if main_source=="agri", mc("34 101 76") ms(D) msize(large)) ///
					(scatter n mean_ss_beta  if main_source=="other", mc("33 126 55") ms(D) msize(large) ///
					ytitle("",axis(1)) yla(none,nogrid axis(1)) ysc(noline axis(1))   ytitle("") xtitle("Correlation between wealth" "and average pollution" "with dust and sea salt removed",size(vlarge)) xlabel(,labsize(large)) ysize(9) xsize(2) legend(off) xline(0))			
					
					
					
					
					graph export "$figures/country_figure_mean_collapsed_edgar_ss.png", as(png) replace

	
clear
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
		
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	